Particle simulator and method of simulating particles

ABSTRACT

The memory size used in DEM calculation of particle having a particle diameter distribution is suppressed. 
     A particle simulator includes a particle-information retaining unit  11  holding particle information including position and velocity information of a particle group, a particle-number changing unit  14  assigning particle numbers specifying particles in an order in accordance with the positions of the particles to the particle, a contact-candidate-list preparing unit  16  selecting particle pairs of a target particle and another particle that may be in contact with the target particle, a contact determining unit  18  calculating contact forces generated between particles in the particle pairs on the basis of particle information and storing the contact forces in the contact-force tables, a contact-force calculating unit  19  extracting contact forces of particles having a diameter greater than particles from contact-force tables using a contact-force reference table  54 , extracting contact forces of particles having a diameter smaller than the particles by specifying the storage positions in the contact-force tables using integrated-contact-candidate numbers s_jgi[i], and calculating the sum of the contact forces, and a particle-information updating unit  20  updating the particle information on the basis of the contact forces of the particle.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to a particle simulator and a method ofsimulating particles.

2. Related Background Art

A discrete element method (DEM), which has originated through researchin civil engineering, has been applied in particle technology. Today,the method is in great demand as a numerical analysis method used tosolve an extremely wide variety of problems associated with not onlyparticles but also continuous subjects that cannot be readily treated.Many other methods of numerically tracking discrete particle motion,such as molecular dynamics (MD), lattice gas automation (LGA), andsmoothed particle hydrodynamics, have also been developed, in additionto the DEM. The DEM has a wide range of application compared to theother numerical particle analysis methods because the calculation ofdetermining interparticle forces in consideration of geometricinformation of the particles, such as the size and shape of theparticles, and solving Newton's motion equation can be performed at highspeed. Size and shape can be considered in methods, such as adiscontinuous deformation analysis (DDA), other than the DEM.Unfortunately, such a method exhibits high initial calculation load, andthus, it is not practical to process a large number of particles.

Commonly known granules and powders (which are hereinafter collectivelyreferred to as “particles”) that can actually be seen, such as sand andflour, usually includes particles of various different sizes unlessspecial processing, such as particle-size adjustment through meshcontrol, is carried out (this state is referred to as a state having a“particle diameter distribution”). Such a particle diameter distributiongreatly influences the homogeneity of the particle and mechanicalproperties, such as deformation and flow. In industries that processparticles, the particle-diameter separation efficiency, theheterogeneity and mixing efficiency of particle layers due to particlesize segregation, and control of transportation and filling is alwaysproblematic. Therefore, to process particles, it is critical toaccurately understand how the effect of a particle diameter distributionon the overall behavior of a particle group.

Since the DEM is a method that can consider the particle size, asdescribed above, it is understandable to use the DEM to determine themicromechanics of particles based on a particle diameter distribution.However, a decrease in the particle diameter to 1/n causes a decrease inthe volume of the particle to (1/n)³. This indicates, n³ particles withreduced volume are required to fill the volume of the original particle(precisely, the number of reduced-volume particles required is less thann³ because there will be gaps between the particles, but the order ofthe number is the same). Therefore, a wide particle diameterdistribution cannot be processed due to restrictions on the amount ofcomputer memory and calculation time. For example, natural soil containsparticles varying in size by approximately five digits, e.g., fromgravels having a diameter of several centimeters to natural clay havinga diameter of several microns. Under an assumption that particles ofeach different size group occupy the same volume, if only one particlehaving a maximum diameter is included, the total number of particleswill be on the order of (10⁵)³.

In the case where a particle diameter distribution is processed usingthe DEM, a disadvantage is imbalance of the number of contact points dueto the difference in particle diameter, in addition to the problem onthe number of particles. For example, if the particle diameter ratio ofthe maximum particle to the minimum particle is r:1, the number ofminimum particle covering the surface of a maximum particle is4(r+1)²/1², whereas the number of maximum particles filling the surfacecovering the surface of a minimum particle is 4(r+1)²/r². In otherwords, the maximum imbalance of the number of contact points is on theorder of r²:1. The imbalance in the number of contact points lead to animbalance in the number of calculation processes performed to determinethe contact forces of the particles and the amount of calculationrequired to determine the sum of the contact forces. Specifically, ifthe particle diameters are not uniform and the particle diameter ratioincreases in the DEM, the number of particles increases simply due togeometric features on the diameter, and causes an imbalance in thenumber of contact points. As a result, the calculation load considerablyincreases. In addition, if a particle diameter distribution exists inthe actual program for the DEM, processing of refined search of contactcandidates of which the contact is to be determined and storingcontact-candidate pairs causes a significantly high calculation loadcompared with a uniform particle size.

In response to such problems concerning a particle diameterdistribution, several methods of DEM programing for increasingprocessing speed associated with extracting and storingcontact-candidate pairs have been proposed (for example, JapaneseUnexamined Patent Application Publication No. 2007-286514 (PatentLiterature 1)).

Today, high speed numerical simulation using GPUs, which have a highcost performance rate compared with CPUs, is actively used. In responseto this, instead of vector processing using expensive super computers,the use of high-speed algorithms for GPUs, which are inexpensive butcapable of super parallel computing equivalent to vector operation, havebeen promoted for particle system simulation (for example, JapaneseUnexamined Patent Application Publication No. 2009-69930 (PatentLiterature 2)).

SUMMARY OF THE INVENTION

The conventional high-speed DEM according to Patent Literature 1involves storing particles in cells and then optimizing the cell sizeand cell search range and is not particularly aimed for a substantialincrease in processing speed by parallelization and vectorization.Specifically, if the process of storing particles into cells by loopingthe particle numbers is performed through single-instructionmultiple-data (SIMD) parallel processing, competitive writing occurs inthe memory when two or more particles are stored in the same cell, whichis an extremely inconvenient situation. Such a problem of competitivewriting in the memory also occurs in particles having a uniform diameterif a plurality of particles is included a single cell. However, if thereis a particle diameter distribution, a distribution of the particlescontained in a single cell is wide in itself, causing an even greaterproblem in increasing the processing speed by parallelization andvectorization.

The detection of contact-candidate particle pairs using a GPU high-speedalgorithm according to Patent Literature 2 is performed by employing amethod of searching cells after storing particles therein in the samemanner as the CPU-based program. However, since competitive writing inthe memory also occurs with the GPU when a plurality of particles isstored in the same cell, several solutions are provided to prevent sucha problem. One solution is a method of repeating several times thestorage of particle and confirmation. Another method use an exclusionfunction (Atomic function) included in the compiler. To preventcompetitive writing in the memory during a routine for calculating theinterparticle contact forces, the same calculation of calculating thecontact forces between particles i and j is performed twice for both theparticles i and the particles j.

However, as described above, if there is a particle diameterdistribution, the number of particles stored in a cell is likely to varygreatly in different cells, and the number of contact particles of alarge particle is significantly larger than that of a small particle.Therefore, the load balance tends to be disturbed. For example, in thecase where GPU calculation is performed for a two-component-systemparticle group using a three-dimensional DEM employing these methods, itis reported that the particle diameter ratio of the two components isset to 1:4, and the calculation time is extended by three times comparedwith a single component system

The DEM simulation of particles considering a particle diameterdistribution leads to a drastic increase in the number of neighboringparticles of which contact-candidate forces are to be considered.Accordingly, there is a problem of increasing a required volume ofmemory to perform DEM calculation.

Unlike in CPU calculation, the memory cannot be increased by hardware inGPU calculation. Accordingly, processing can only be performed using alimited volume of memory in a graphic card. The shortage in GPU memorymay be supported by the CPU memory through transfer communication.However, this is not recommended since the transfer speed between theGPU memory and the CPU memory is extremely low, and the high calculationefficiency of the GPU is not fully used. Since the memory size for anexisting GPU used in DEM calculation is 4 GB at a maximum, the problemconcerning the memory size used prevents the use of GPU in particlesimulation. This problem exists not only in GPUs but also in hardware(for example, a vector processor) (hereinafter referred to as“shared-memory parallel computer”) that performs shared-memory SIMDparallel processing similar to that performed by a GPU.

The present invention has been conceived to solve the problem describedabove and provides a particle simulator and a method of simulatingparticles that are capable of reducing the amount of memory during DEMcalculation of particle having a particle diameter distribution by ashared-memory parallel computer.

To solve the problem described above, a particle simulator according tothe present invention simulating the behavior of particles in a particlegroup having a diameter distribution in a work area through calculationof positions and velocities of the particles based on values of contactforces generated between the particles and other particles in theparticle group, the particle simulator comprising: input means forinputting particle information including position information andvelocity information of the particles; particle-information storingmeans for storing the particle information; assigning means forassigning specific information items to the particles to identify theparticles in an order in accordance with the position information;contact-candidate selecting means for selecting target particles fromparticle group in the ascending order of the specific information itemsand selecting non-target particles paired with the target particles ascontact-candidate pairs, the non-target particles presumably being incontact with the target particles and having a diameter smaller than thediameter of the target particles; contact-force calculating means forcalculating the values of contact forces generated between the twoparticles in the contact-candidate pairs selected by thecontact-candidate selecting means based on the position information andthe velocity information of the paired particles and storing informationitems of the calculated values of contact forces in an ordercorresponding to the order of the contact-candidate pairs incontact-force tables; contact-force-reference-table preparing means forpreparing contact-force-reference tables including reference informationused to refer to the contact-force tables to retrieve the values ofcontact forces generated between the particles and neighboring particlesneighboring the particles and having a diameter greater than thediameter of the particles; integral-contact-candidate-number calculatingmeans for calculating integral contact-candidate numbers being integralvalues of contact-candidate numbers indicating the numbers ofneighboring particles neighboring the particles having a diametersmaller than the diameter of the particles; contact-force-sum computingmeans for extracting the values of contact forces generated between theparticles and other particles having a diameter greater than thediameter of the particles by using the reference information in thecontact-force-reference tables, extracting the values of contact forcesgenerated between the particles and other particles having a diametersmaller than the diameter of the particles by specifying the storagepositions of the values of contact forces in the contact-force tablesusing the integral contact-candidate numbers, and computing the sum ofcontact forces for each particles using the extracted values of contactforce; calculating means for calculating new position information andvelocity information of the particles based on the values of contactforces of the particles calculated by the contact-force-sum computingmeans; and particle-information updating means for updating the particleinformation stored in the particle-information storing means using thenew position information and velocity information calculated by thecalculating means.

Similarly, to solve the problem described above, a method of simulatingparticles in a particle simulator of the invention simulating thebehavior of particles in a particle group having a diameter distributionin a work area by calculating the positions and velocities of theparticles based on values of contact forces generated between theparticles and other particles in the particle group, the methodcomprising, the steps of inputting particle information includingposition information and velocity information of the particles andstoring the particle information in particle-information storing means;assigning specific information items to the particles to identify theparticles in an order in accordance with the position information;selecting target particles from particle group in the ascending order ofthe specific information items and selecting non-target particles pairedwith the target particles as contact-candidate pairs, the non-targetparticles presumably being in contact with the target particles andhaving a diameter smaller than the diameter of the target particles;calculating the values of contact forces generated between the twoparticles in the contact-candidate pairs selected in the selecting stepbased on the position information and the velocity information of thepaired particles and storing information items of the calculated valuesof contact forces in an order corresponding to the order of thecontact-candidate pairs in contact-force tables; preparingcontact-force-reference tables including reference information used torefer to the contact-force tables to, retrieve the values of contactforces generated between the particles and neighboring particlesneighboring the particles and having a diameter greater than thediameter of the particles; calculating integral contact-candidatenumbers being integral values of contact-candidate numbers indicatingthe numbers of neighboring particles neighboring the particles having adiameter smaller than the diameter of the particles; extracting thevalues of contact forces generated between the particles and otherparticles having a diameter greater than the diameter of the particlesby using the reference information in the contact-force-referencetables, extracting the values of contact forces generated between theparticles and other particles having a diameter smaller than thediameter of the particles by specifying the storage positions of thevalues of contact forces in the contact-force tables using the integralcontact-candidate numbers, and computing the sum of contact forces foreach particles using the extracted values of contact force; calculatingnew position information and velocity information of the particles basedon the values of contact forces of the particles calculated in theextracting step; and updating the particle information stored in theparticle-information storing means using the new position informationand velocity information calculated by the calculating step.

According to such a particle simulator and a method of simulatingparticles, for a particle neighboring a particle and having a diametersmaller than that of the particle among all particles in the work area,only the integral contact-candidate numbers, which represent the numbersof particles, is stored, and the contact forces are extracted byspecifying the storage positions of the contact-force tables using theintegral contact-candidate numbers. In this way, the contact-forcereference tables do not have to contain information associated withparticles having a diameter smaller than that of particles amongparticles neighboring the particles, and merely information of particleshaving a diameter greater than that of the particles needs to be stored.During simulation on a particle group having a particle diameterdistribution, if particles having a diameter smaller than that of theparticle are considered, the number of particles that may simultaneouslycontact the particle will increases. Thus, a large memory area must beprovided in advance for the contact-force reference table. In contrast,the configuration according to the present invention requiresconsideration of only the contact between the particle and particleshaving a diameter greater than that of the particle. Thus, the number ofparticles for which contact determination need to be performed isdecreased, and the memory region required for preparing thecontact-force reference tables can be decreased. As a result, the memorysize used can be reduced. Such a reduction in the memory size usedenables particle simulation using a shared-memory parallel computer,such as a GPU.

In the particle simulator, it is desirable that the work area besegmented into a plurality of cells, each cell having a cell number, theparticle information include the cell numbers of the cells specified bythe position information and containing the particles, and the assigningmeans identifies the particles in an order corresponding to the cellnumbers and assigns specific information items specified in accordancewith the order corresponding to the particle diameter to particles inthe same cell.

Similarly, in the method of simulating particles it is desirable thatthe work area be segmented into a plurality of cells, each cell having acell number; the particle information includes the cell numbers of thecells specified by the position information and containing theparticles; and in the assigning step, the particles in an ordercorresponding to the cell numbers are identified and assigns specificinformation items in accordance with the order corresponding to theparticle diameter to particles in the same cell.

According to this configuration, since specific information is assignedby identifying the particles in the order of the cell number, andidentifying the particles in the same cell in an order in accordancewith the particle diameter, the particles are sorted by cell number, andparticles in the same cell are sorted by particle diameter. Uponpreparation of the contact-force reference tables, particles having asmall diameter can be easily recognized on the basis of the specificinformation and do not have to be selected for entry to thecontact-force reference tables. As a result, the calculation efficiencycan be improved during preparation of the contact-force reference table.

In the particle simulator, it is desirable that the contact-candidateselecting means selects the target particles from the particle group inthe ascending order of the specific information items and preparescontact-candidate lists of the target particles including, for all ofthe particles in the work area, contact-candidate-pair informationincluding pairs of the target particles and the non-target particlespresumably being in contact with the target particles and having adiameter smaller than the diameter of the target particles andcontact-candidate numbers assigned to the pairs in the ascending order;the contact-force calculating means calculates the values of contactforces generated between the paired particles corresponding to thecontact-candidate-pair information included in the contact-candidatelists are calculated based on the position information and the velocityinformation of the paired particles and store contact-force informationlinking the calculated values of contact forces and thecontact-candidate number in contact-force tables; thecontact-force-reference-table preparing means prepares thecontact-force-reference tables linking the contact-candidate numbersassociated with the particles and other particles with the particles andthe other particles as reference information used to refer to thecontact-force tables to retrieve the values of contact forces generatedbetween the particles and the other particles presumably being incontact with the particles and having a diameter greater than thediameter of the particles; the integral-contact-candidate-numbercalculating means calculates integral contact-candidate numbers s_jgi[i](where i represents specific information items of the particles) beingintegral values integrating, in the order of the specific informationitems, contact-candidate numbers indicating the number of otherparticles presumably being in contact with the particles and having adiameter smaller than the diameter of the particles; and thecontact-force-sum computing means extracts the values of contact forcesof the particles from the contact-force tables using the referenceinformation of the contact-force-reference tables, extract the values ofcontact forces corresponding to the contact-candidate numbers s_jgi[i−1]to s_jgi[i]−1 from the contact-force tables using the integralcontact-candidate numbers s_jgi[i], and calculate the sum of the contactforces of the particles using the extracted information of the contactforces.

The particle simulator and the method of simulating particles accordingto the present invention, can reduces the memory size used during DEMcalculation of particles having a particle diameter distribution by ashared-memory parallel computer.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a functional block diagram of a particle simulator accordingto an embodiment of the present invention.

FIG. 2 is a schematic view of a Voigt model applied to a particle modelaccording to this embodiment.

FIG. 3 illustrates a work area according to this embodiment.

FIG. 4 illustrates a flow of changing particle numbers.

FIG. 5 illustrates the process of changing particle numbers and storingthe maximum value and minimum value of the particle number in each cell.

FIG. 6 illustrates an example configuration of a neighboring-particletable.

FIG. 7 illustrates an example configuration of a contact-candidate list.

FIG. 8 illustrates a process of confirming whether the paired particlesin the current contact-candidate list were paired in the previous step.

FIG. 9 illustrates an example configuration of a contact-force referencetable.

FIG. 10 illustrates the correspondence among a contact-force referencetable, integral-contact-candidate numbers, and contact-force boxes.

FIG. 11 is a flow chart illustrating the process performed by theparticle simulator according to this embodiment.

FIG. 12 is a diagram illustrating a comparison of the memory size usedin this embodiment and a comparative example, when the ratio of themaximum particle diameter to the minimum particle diameter is varied inthe particle model.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

Embodiments of the present invention will be described below withreference to the accompanying drawings. The same components arerepresented by the same reference numerals without repeated description.

FIG. 1 is a functional block diagram of a particle simulator 10according to an embodiment of the present invention. As illustrated inFIG. 1, the particle simulator 10 includes a particle-informationretaining unit (input means and particle-information storing means) 11,a particle-information acquiring unit 12, acontact-candidate-list-update determining unit 13, a particle-numberchanging unit (assigning means) 14, a neighboring-particle-tablepreparing unit 15, a contact-candidate-list preparing unit(contact-candidate-particle selecting means andintegrated-contact-candidate-number calculating means) 16, acontact-force-reference-table preparing unit(contact-force-reference-table preparing means) 17, a contactdetermining unit (contact-force calculating means) 18, a contact-forcecalculating unit (contact-force-sum computing means) 19, and aparticle-information updating unit (calculating means andparticle-information updating means) 20

The particle simulator 10 is a computer including hardware such as agraphics processing unit (GPU), a random access memory (RAM) and a readonly memory (ROM), which are main storage devices, an auxiliary storagedevice, such as a hard disk drive, an input device, such as an inputkeyboard, an output device, such as a display, and a communicationmodule, which is a data transmission/reception device, such as a networkcard. The functions of the particle simulator 10 described below withreference to FIG. 1 are achieved by predetermined computer softwareloaded into the hardware, such as the GPU and the RAM to operate theinput device, the output device, and the communication module under thecontrol of the GPU and to read/write data from/to the RAM and theauxiliary storage device.

A particle model used in this embodiment will be described below. Theparticle simulator 10 according to this embodiment performs a DEMsimulation of a particle group having a certain particle-diameterdistribution. More specifically, a particle model including particleshaving various diameters is used.

In such a particle model, the normal component F_(n) and tangentialcomponent F_(t) of an interparticle contact force generated betweenindividual particles are represented by the following expressions byusing a Voigt model 71, which is illustrated in FIG. 2.

F _(n)=η_(n) v _(n) +K _(n) x _(n)  (1)

F _(t)=min[η_(t) v _(t) +K _(t) x _(t),μ_(t) |F _(n) |n _(t)]  (2)

where x represents a relative displacement vector, v represents arelative velocity vector, K represents the elastic modulus of a linearspring 72, and η represents the viscous damping coefficient of a dashpot73. The subscripts n and t represent a normal component and a tangentialcomponent, respectively. μ_(t) represents the coefficient of slidingfriction, and min represents a function representing the value of thesmaller of input values. The relationship between the viscous dampingcoefficient η and a coefficient of restitution e is represented by thefollowing expression.

$\begin{matrix}{\eta_{q} = {{- \frac{2\ln \; e}{\sqrt{\pi^{2} + {\ln^{2}e}}}}\sqrt{\frac{2m_{i}m_{j}}{m_{i} + m_{j}}K_{q}}}} & (3)\end{matrix}$

where m_(i) and m_(j) represent the masses of particles i and j,respectively. The subscript q is equals to n for a normal component andt for a tangential component.

Accordingly, equations of motion of translation and rotation around thecenter of gravity, which are processed in the DEM, are represented belowby Expressions 4 and 5.

$\begin{matrix}{{m_{p}\frac{\partial u_{p}}{\partial t}} = {F_{n} + F_{t} - {m_{p}g}}} & (4) \\{{I_{p}\frac{\partial\omega_{p}}{\partial t}} = {{R \times F_{t}} - {\mu_{r}b{F_{n}}n_{r}}}} & (5)\end{matrix}$

where u_(p) represents the velocity of a particle, ω_(p) represents theangular velocity of the particle, m_(p) represents the mass of theparticle, I_(p) represents the moment of inertia of the particle, and Rrepresents a position vector from the center of the particle to acontact point. The subscript p is equals to i for a particle i and j fora particle j. The second term on the right side indicates the rollingfriction, where μ_(r) represents the coefficient of rolling friction andb represents the width of the contact surface. Expressions 1 and 2 ofthe Voigt model 71 are used by the contact determining unit 18 tocalculate the interparticle contact force, whereas the equations ofmotion, i.e., Expressions 4 and 5, are used by the particle-informationupdating unit 20 to calculate the coordinates of a particle.

In this embodiment, a work area 51, which is the area in which theparticles move, is divided into cubic cells having sides, each having alength of csize, as illustrated in FIG. 3. The cell size csize isdetermined by the following expression using a parameter φ(0≦φ≦1).

csize=2.0{(1.0−φ)r _(max) +φr _(min))}(1.0+α)  (6)

where r_(max) represents the maximum particle diameter in the particlemodel arranged in the work area, r_(min) represents the minimum particlediameter, and α represents a parameter for adjusting the updatingfrequency, which is set to, for example, 0.2. Each cell in the work area51 is assigned a cell number cell_id. The cell number isone-dimensional, as illustrated in FIG. 3, and is represented by, forexample, the fowling expressions.

cell_(—) id=i.x+id.x×i.y+id.x×id.y×i.z

i.m={0, 1, 2, . . . , id.m−1}(m=x, y, z)  (7)

where cell_id represents a cell number, i.m represents the cell positionin the m direction, id.m represents the number of cells in the mdirection. After simulation conditions, i.e., the maximum particlediameter r_(max), the minimum particle diameter r_(min), and theparameters φ and α, are set, the particle simulator 10 automaticallycalculates and sets the cell size csize of the work area 51 and assignscell numbers to the defined cells.

The particle-information retaining unit 11 holds particle information ofeach particle in the work area 51. The particle information includescoordinates (position information), translation velocity (velocityinformation), rotational velocity (velocity information), particleradius, mass, particle number, and cell number. The coordinates arethree-dimensional coordinates (pos.x, pos.y, pos.z) in thethree-dimensional work area 51, which is illustrated in FIG. 3. Theparticle number is the number for identifying a particle. The cellnumber is the cell number assigned to the cell containing the particle.To begin simulation, the particle-information retaining unit 11 sets theparticle radii and masses of particles in accordance with theparticle-diameter distribution used in the simulation and randomlyassigns the coordinates, velocities, particle numbers to theseparticles. Each cell number corresponding to each particle is assignedthe cell number of the cell contain the particle that is identified onthe basis of the coordinates randomly assigned to the particle. Eachinformation item generated in this way is input and stored in theparticle-information retaining unit 11. When the particle informationitems are updated by the particle-information updating unit 20, which isdescribed below, during the particle simulation, theparticle-information retaining unit 11 replaces the stored particleinformation items with the particle information items updated by theparticle-information updating unit 20.

The particle-information acquiring unit 12 acquires particle informationfrom the particle-information retaining unit 11. Theparticle-information acquiring unit 12 acquires initial particleinformation set in advance from the particle-information retaining unit11 at the beginning of the particle simulation. Every time theparticle-information retaining unit 11 updates the particle informationduring the particle simulation, the particle-information acquiring unit12 acquires the updated particle information from theparticle-information retaining unit 11. The particle-informationacquiring unit 12 sends the acquired particle information to thecontact-candidate-list-update determining unit 13.

The contact-candidate-list-update determining unit 13 determines whetheror not to update a contact-candidate list 53, which is described below.The contact-candidate list 53 is updated, for example, if a particle ofwhich the integrated traveling distance is greater than rα is present,where r represents the particle radius and α represents a parameter foradjusting the updating frequency, which is set to, for example, 0.2. Theintegrated traveling distance is determined for each particle, forexample, every time the particle-information acquiring unit 12 acquiresupdated particle information from the particle-information retainingunit 11. The integrated traveling distance of a particle is determinedby integrating the latest travelling distance of the particle, which iscalculated from the difference between coordinate data in the particleinformation acquired in the current step and the coordinate data in theparticle information acquired in the last step. Thecontact-candidate-list-update determining unit 13 holds the particleinformation acquired in the last step, the parameter α, and theintegrated traveling distance of particles. Upon reception of updatedparticle information from the particle-information acquiring unit 12,the contact-candidate-list-update determining unit 13 updates theintegrated traveling distance and determines the existence of a particleof which the updated integrated traveling distance is greater than rα.If a particle of which the integrated traveling distance is greater thana predetermined value is present, the contact-candidate-list-updatedetermining unit 13 determines that the contact-candidate list should beupdated and transmits the updated particle information to theparticle-number changing unit 14. The integral value of the particletraveling distances is initialized every time the contact-candidate list53 is updated. If the contact-candidate-list-update determining unit 13determines that the contact-candidate list does not have to be updatedbecause there are no particles of which the integrated travelingdistance is greater than rα, the contact-candidate-list-updatedetermining unit 13 does not send particle information to theparticle-number changing unit 14 and sends the updated particleinformation to the contact determining unit 18.

The particle-number changing unit 14 identifies the particles in theorder of the cell numbers (i.e., in the order corresponding to theposition information) and assigns particle information (specifyinginformation) to particles contained in the same cell in order of theparticle diameter for each of the particles in the work area 51. Inother words, the particle-number changing unit 14 renumbers theparticles on the basis of the cells containing the particles in the workarea 51 and the particle diameter and sorts out the particles by the newparticle numbers.

Upon receiving the particle information from thecontact-candidate-list-update determining unit 13, the particle-numberchanging unit 14 creates information sets, each linking the particlenumber of individual particles with the diameter of the particle and thecell number of the cell containing the particle, as illustrated in FIG.4, and stores these information sets in the order of the particlenumbers (FIG. 4( a)). The information sets are rearranged in theascending order of cell numbers (4(b)), and then the particle numbersare reassigned in the order of cell numbers (FIG. 4( c)). The particlediameters of particles having the same cell number are compared torearrange these particles in the descending order of the particlediameter (FIG. 4( d)). Then, the particle numbers are reassigned on thebasis of particle diameter (FIG. 4( e)).

Through this rearrangement process of the particle numbers, for example,the particle numbers of particles disposed in the work area 51, asillustrated in FIG. 5( a), are renumbered according to the cell numbersand particle diameters, as illustrated in FIG. 5( b). Theparticle-number changing unit 14 transmits the particle informationcontaining the change in particle numbers to theneighboring-particle-table preparing unit 15.

The neighboring-particle-table preparing unit 15 preparesneighboring-particle tables 52 containing information about particlesneighboring individual particles in the work area 51. In theneighboring-particle-table preparing unit 15, “neighboring particles”are defined as particles in adjoining cells that adjoin a target cell,which contains a target particle, and in neighboring cells, whichoutwardly adjoin the adjoining cells, having a particle diameter greaterthan that of the target particle or having a particle number, which isassigned in the initial position, greater than that of the targetparticle in the case where the particle diameters are the same. Theassignment of neighboring particles will be described in detail below.

Memory regions are reserved in advance in the neighboring-particle-tablepreparing unit 15 for preparing the neighboring-particle tables 52. Thememory regions for neighboring-particle tables are set to the same sizein consideration of the number of neighboring particles that may comeinto contact with each particle.

Upon receiving the particle information from the particle-numberchanging unit 14, the neighboring-particle-table preparing unit 15stores the minimum particle number Pnum_in_cellmin and the maximumparticle number Pnum_in_cellmax selected among particles having the samecell number (FIG. 5( c)). In FIG. 5( c), the maximum particle number isstored in the upper left region and the minimum particle number isstored in the lower right region in a cell. More specifically, each cellcell_id is provided with a sequence in which information about theminimum particle number Pnum_in_cellmin[cell_id] and the maximumparticle number Pnum_in_cellmax[cell_id] of particles contained in thecell is stored. In this way, each cell of the cell number cell_idcontains particles numbered from Pnum_in_cellmin[cell_id] toPnum_in_cellmax[cell_id] in the descending order of particle diameter.

The cell number nei_cell_id of a cell adjoining the cell cell_id isrepresented by the following expression.

nei_cell_(—) id[n]=cell_(—) id+con[n]  (8)

where con represents a constant determined by id.x and id.y. Since theparticle numbers are sorted by the cell number, each neighboring cellcontains particles continuously numbered fromPnum_in_cellmin[nei_cell_id] to Pnum_in_cellmax[nei_cell_id].

The neighboring-particle-table preparing unit 15 prepares storage boxes(sequence) nei_pn[i][box_id] for storing neighboring particles of eachparticle, such as that illustrated in FIG. 6. Among the particles in asearch range (which is defined below) around a particle i, particlesthat satisfy the following Conditions a, b, and c are stored in theboxes in the ascending order of the box numbers (box_id): a) particlesof which the interparticle distance to the particle i is smaller thandis_μm; b) particles having a diameter greater than the diameter of theparticle i; and c) if particles have the same diameter as the diameterof the particle i, the particles having a particle number (i=p_id)assigned in the initial position greater than the particle number of theparticle i, where dis_lim=(r_(i)+r_(j))(1.0+α), r_(i) and r_(j)represent the particle radii of paired neighboring particles i and j,respectively, and α represents a parameter for adjusting the updatingfrequency of contact candidates. For example, α is set to 0.2. Theneighboring-particle-table preparing unit 15 calculates dis_lim usingthe particle radii r_(i) and r_(j), which are included in the particleinformation received from the particle-number changing unit 14, and theparameter α which is stored in advance in the neighboring-particle-tablepreparing unit 15.

A particle that satisfies Conditions (a), (b), and (c) is referred to asa “neighboring particle” of particle i. The collection of the series ofboxes, which is illustrated in FIG. 6, provided for every particle i inthe work area 51 is referred to as a “neighboring-particle table 52” inthis embodiment. The neighboring-particle table 52 is a two-dimensionalsequence nei_pn[i][box_id] associated with each particle i and thenumber (box_id) of neighboring particles of each particle i.

Through this procedure, for target particles i, only particles j havinga particle diameter greater than that of the particle i are stored inthe neighboring-particle table 52. Therefore, increasing the memory arearequired for the neighboring-particle table 52 is inhibited.

In order to determine the neighboring particles of the particle i, theneighboring-particle-table preparing unit 15 has to search for particlesat a distance (r_(max)+r_(i))(1.0+α) from particle i. In addition to the26 cells adjoining the cell containing the particle i, the search rangeis expanded to outer cells in accordance with the radius r_(max) of thelargest particle in the particle model.

The range of neighboring cells searched by theneighboring-particle-table preparing unit 15 is determined using thefollowing expression, where csize represents the cell size calculated byExpression 6, r_(max) represents the maximum particle radius of thelargest particle in the work area, and r_(i) represents the radius ofthe target particle i.

nei_cell=2·ceil{2.0(r _(max) +r _(i))(1.0+α)/csize}−1  (9)

where ceil represents a function for rounding up the fractional part andis stored in advance in the neighboring-particle-table preparing unit15, and nei_cell represents the number of cells to be searched in theplus and minus directions of the x, y, and z directions from the cellcontaining the particle i.

Consequently, as the cell size is reduced, there is a worry that thenumber of cells to be searched increases exponentially, increasing thecalculation load. Meanwhile, for the cells on the outer side of theadjoining 26 cells, the order of the particles sorted in accordance withparticle diameter can be used to increase the efficiency of searchingneighboring particles. Specifically, each cell contains particles fromcell numbers Pnum_in_cellmin[cell_id] to Pnum_in_cellmax[cell_id] in thedescending order of particle diameter, and determination of whether ornot a particle is a neighboring particle is also performed on particlesin the descending order of particle diameter. Since the diameter ofparticles that do not satisfy above-described Condition (a) isdetermined on the basis of the cell containing the particle i and thepositional relationship between the particle i and the neighboringparticles, the interparticle distance does not have to be checked inevery cell up to Pnum_in_cellmax[cell_id], but the search can beterminated at a particle number corresponding to a particle having adiameter that does not satisfy Condition (a). A particle diameter r_(j),which does not satisfy Condition (a) for the particle i, satisfies thefollowing expression:

2ceil{2.0(r _(j) +r _(i))(1.0+α)/csize}−1<dis cell  (10)

where, dis cell represents the number of cells from the cell containingthe particle i to the cell containing particle j in one of the x, y, andz directions having the least cell number.

The neighboring-particle-table preparing unit 15 stores neighboringparticle numbers n_jli[i] of particles neighboring the particle i andparticle numbers and n_jgi[i] of particles that satisfy Condition (a)but do not satisfy Conditions (b) and (c) in correspondingone-dimensional sequences.

In FIG. 6, search is performed for a particle numbered 4 (i=4). If theparticles numbered 5, 6, 7, and 8 satisfy Condition (a) as the searchresult, if the particles numbered 7 and 8 are neighboring particles, andif the particles numbered 5 and 6 do not satisfy Conditions (b) and (c),the following equations hold:

nei _(—) pn[4][0]=7

nei _(—) pn[4][1]=8

n _(—) jgi[4]=2

n _(—) jli[4]=2

When the neighboring-particle-table preparing unit 15 receives particleinformation from the particle-number changing unit 14, cells in thesearch range nei_cell are defined for every target particle i, and theparticle numbers of particles other than the particle i in the searchrange that satisfy Conditions (a), (b), and (c) are stored in thesequence nei_pn[i][box_id] of the neighboring-particle table 52. Thenumber of neighboring particles is stored in sequence n_jli[i], whereasthe number of particles that satisfy Condition (a) but do not satisfyConditions (b) and (c), i.e., the number of particles in the searchrange and not being stored in the neighboring-particle table 52 becausetheir diameters are smaller than the diameter of the target particle iis stored in sequence n_jgi[i]. The information of nei_pn[i][box_id],n_jli[i], and n_jgi[i] is sent together with particle information to thecontact-candidate-list preparing unit 16.

In this embodiment, only information of particles that have a diametergreater than or equal to the diameter of the target particle is storedin the neighboring-particle table 52. In a multi-component systemincluding particles having various different diameters, the number ofsmall particles that are contact candidates for large particles may beenormous. That is, numerous numbers of the particles having a diametersmaller than that of the target particle can simultaneously contact thetarget particle. Therefore, if particle having a diameter smaller thanthat of the target particle are included in the neighboring particles,number of sequences for the neighboring-particle table 52 has to beincreased in advance in consideration of a maximum. In such a case, anenormous amount of memory for the sequences in the neighboring-particletable 52 is required, and the global memory of the GPU may be exhausted.In this embodiment, since only particles that are larger than the targetparticle are stored in the neighboring-particle table 52, the amount ofmemory used for the sequences prepared in advance for theneighboring-particle table 52 is inhibited.

The contact-candidate-list preparing unit 16 sequentially selects targetparticles in the ascending order of the particle number (specificinformation) from among all particles in the work area 51 and selectscontact candidate pairs of target particles and particles that may be incontact with the target particles and that have a diameter smaller thanthat of the target particles. Specifically, the contact-candidate-listpreparing unit 16 prepares the contact-candidate list 53 containinginformation on particle pairs of particles that are selected from allparticles in the work area 51 and may be in contact with each other.More specifically, the contact-candidate list 53 includesone-dimensional sequences list_i[Lnum] and list_j[Lnum] in which theparticle numbers of particles i and j, which may be in contact with eachother, are respectively stored, where Lnum represents acontact-candidate number given to each pair of particles in contact witheach other.

Upon receiving the information (neighboring-particle table 52nei_pn[i][box_id], neighboring particle numbers n_jli[i] and n_jgi[i],and particle information) from the neighboring-particle-table preparingunit 15, the contact-candidate-list preparing unit 16 determines aprefixed sum s_jgi of n_jgi, as illustrated in FIG. 7, using thefollowing expression.

$\begin{matrix}{{{s\_ jgi}\lbrack i\rbrack} = {\sum\limits_{{p\_ id} = 0}^{i}{{n\_ jgi}\lbrack{p\_ id}\rbrack}}} & (11)\end{matrix}$

The prefixed sum s_jgi[i] is a one-dimensional sequence representing anintegral value, which integrates, up to the particle number i, theparticle numbers of particles not stored in the neighboring-particletable 52 of particle i, i.e., the number of particles n_jgi[i] (numberof contact candidates) having a diameter smaller than that of theparticle i among all particles existing in a predetermined search rangearound the particle i (a neighboring cell defined by Expression 9). Inthis embodiment, the prefixed sum s_jgi[i] is referred to as “integralcontact-candidate number.” Specifically, the integral contact-candidatenumber s_jgi[i] is a value integrating the contact candidates n_jgi[i]in the order of particle numbers.

The contact-candidate-list preparing unit 16 sequentially selects targetparticles i in the ascending order of the particle numbers, re-searchesthe cell, which is the searching range around the particle i defined byExpression 9, and extracts the corresponding particles j, which satisfythe following conditions, from among particles contained in thesearching range around the particle i: a) particles of which theinterparticle distance to the particle i is smaller than dis_lim; b)particles having a diameter smaller than the diameter of the particle i;and c) if particles have the same diameter as the diameter of theparticle i, the particles having a particle number assigned in theinitial position smaller than the particle number of the particle i.Contact-candidate numbers Lnum are assigned to the pairs of particles iand j in the ascending order, and the particle number of the particles iand j are respectively stored in the one-dimensional sequences list_iand list_j to prepare the contact-candidate list 53.

Satisfying Conditions a), b), and c) is equivalent to satisfyingabove-described Conditions (a) of the extracting conditions for theneighboring particles to be stored in the neighboring-particle table 52but not satisfying Conditions (b) and (c). The particle pairs not storedin the neighboring-particle table 52 are the contact-candidate pairs(for example, the contact-candidate pairs are pairs including particles5 and 6 (i=5 and i=6), if the target particle i is numbered 4 (i=4)).

For the example illustrated in FIG. 7, the content of thecontact-candidate list 53 is prepared as Table 1.

TABLE 1 Lnum 0 1 2 3 4 5 6 7 list_i[Lnum] 1 4 4 5 7 7 7 8 list_j[Lnum] 25 6 6 4 5 6 4

Specifically, the contact-candidate list 53 containscontact-candidate-pair information including pairs of particles list_iand list_j, which may be in contact with each other for all particles inthe work area 51, together with particle numbers Lnum.

The amount of expansion/contraction of the spring required fordetermination of the contact force in the normal direction is anintegral value. If paired contact-candidate particles are in contactwith each other both in the previous step and the current step, theamount of expansion/contraction of the spring in the normal-directionshould be retained for updating the contact-candidate pair list. Thepaired contact candidates Blist_i and Blist_j in the previous step arethus searched to determine whether the paired contact-candidates list_iand list_j in the current step were also pairs in the previous step. Ifthey were pairs, the normal vector and the spring force (the amount ofexpansion/contraction of the spring) in the tangential direction in theparticle model illustrated in FIG. 2 are stored in sequences of newcandidate list numbers.

Parameters and sequences in the previous step and the current step arerepresented as listed in the table below.

Previous step Current step Particle number Bi = Bpn[i], Bj = Bpn[j] i, jContact-candidate list Blist_i, Blist_j list_i, list_j Integralcontact-candidate Bs_jgi s_jgi number Normal vector and contact BSpringSpring force in tangential direction

If a contact-candidate pair of the particle numbers i=list_i[Lnum] andj=list_j[Lnum] is present, the particle numbers in the previous step isBi=Bpn[i] and Bj=Bpn[j], respectively. The contact candidatesBlist_j[BLnum] within the range of Bs_jgi[Bi−1] to Bs_jgi[Bi]−1 arechecked to confirm whether they match Bj. If a contact candidate matchesBj, the corresponding list number BLnum is used to substitute the normalvector and contact force BSpring[BLnum] in the tangential direction inthe previous step into the current list-number sequence Spring[Lnum].

A process of determining whether paired particles 4 and 6 (i=4, j=6) inthe current contact-candidate list numbered Lnum=2 were also a pair inthe previous step will be described with reference to FIG. 8.

The particle numbers of particles i=4 and j=6 in the previous step areBi=2=Bpn[4] and Bj=3=Bpn[6], respectively, which are determined usingthe Bpn[ ] function. Similar to the relationship between i and j, therelationship between Bi and Bj always satisfy Conditions a), b), and c)of the contact-candidate list 53 for the following reason. Conditionsa), b), and c) are based on the particle diameter and particle number inthe initial position, which do not change by sorting. Accordingly, theparticle numbers of contact candidates of Bi are stored inBlist_j[BLnum] having list numbers BLnum within the range ofBs_jgi[Bi−1] and Bs_jgi[Bi]−1 (which, in this example, is a range from 1to 1). Blist_j[BLnum] is checked against Bj within the range of BLnum.When BLnum=1, Blist_j [1]=3, which is also consistent with Bj=3.

Accordingly, the candidate pair in the current list numbered Lnum=2 wasa pair in the previous step, having a list number BLnum=1. The contactcondition in the previous step can be retained by substituting thenormal vector and the contact force in the tangential direction ofBLnum=1 into the sequence Lnum=2 (Spring[2]=BSpring[1]).

Upon receiving various items of information (neighboring-particle table52 nei_pn[i][box_id], numbers of neighboring particles n_jli[i] andn_jgi[i], and particle information) from the neighboring-particle-tablepreparing unit 15, the contact-candidate-list preparing unit 16 assignsa number Lnum to each pair of a target particle i and a particle j,which satisfies above-described Conditions a), b), and c), among theparticles in the searching range and stores the particle numbers of thepaired particles i and j in the sequences list_i[Lnum] and list_j[Lnum], respectively, in the contact-candidate list 53. Furthermore, anintegral contact-candidate number s_jgi[i], which is an integral valueof the numbers n_jgi[i] of particles present in the cells surrounding aparticle i and having a diameter smaller than that of the particle i iscalculated for each target particle i. Pairs, which were contactcandidate pairs in the previous step are extracted from the currentcontact-candidate pairs list_i[Lnum] and list_j [Lnum], and the contactforce of these extracted pairs is stored in Spring[Lnum] by retainingthe values from BSpring[BLnum]. Then, the derived contact-candidate list53 list_i[Lnum] and list_j [Lnum] and the neighboring-particle table 52nei_pn[i][box_id] are sent to the contact-force-reference-tablepreparing unit 17, and the information of the contact-candidate list 53list_i[Lnum] and list_j[Lnum], the integral contact-candidate numbers_jgi[i], and the retained contact force Spring[Lnum] of thecontact-candidate pairs, which were also pairs in the previous step, aresent, together with the particle information, to the contact determiningunit 18.

The contact-force-reference-table preparing unit 17 generates a functionRef[i] [box_id] to calculate candidate list numbers Lnum fromcorresponding box numbers box_id. In this embodiment, a two-dimensionalsequence Ref[i][box_id] representing this function is referred to as“contact-force reference table 54.” As illustrated in FIG. 9, thecontact-force reference table 54 collects references (Lnum) (referenceinformation) for every particle i (i=4 in FIG. 9) on the contact forceFij[Lnum] of the particle i in the contact-candidate list 53 thatreceives from particle j. Specifically, in this embodiment, thecontact-force reference table 54 includes only references (Lnum) forcontact forces between particles i and neighboring particles stored inthe neighboring-particle table 52 associated with the target particlesi, i.e., particles having a diameter greater than that of the targetparticles i. The references for contact forces (Lnum) arecontact-candidate-pair-list numbers (contact-candidate numbers) Lnumused for obtaining the respective contact forces from contact-forceboxes each containing two different one-dimensional sequences Fijt[Lnum]and Fijr[Lnum] respectively holding the translational component androtational component of an interparticle contact force, which isdescribed below, and a contact force Fij[Lnum] including β[Lnum], whichis used in Expression 13 described below.

Memory regions are reserved in advance in thecontact-force-reference-table preparing unit 17 for preparing thecontact-force reference table 54. The memory regions for thecontact-force reference table 54 are set to the same size inconsideration of the number of neighboring particles that may come intocontact with each particle.

Upon receiving the contact-candidate list 53 list_i[Lnum] andlist_j[Lnum] and the neighboring-particle table 52 nei_pn[i][box_id]from the contact-candidate-list preparing unit 16, thecontact-force-reference-table preparing unit 17 searches theneighboring-particle table 52 of particles j for a particle i(i=list_i[Lnum]) and a particle j (j=list_j[Lnum]), which are the Lnumthparticle pair in the contact-candidate list 53, to acquire theneighboring-particle-table number k (k=box_id) containing the particlei. As described above, since a particle j of a contact-candidate pairhas a diameter smaller than that of the particle i, the particle i isstored in the neighboring-particle table 52 of the particle j as aneighboring particle. The contact-candidate-list number Lnum of theparticles i and j is stored in the kth sequence Ref[j][k] of thecontact-force reference table 54 of the particle j.

The contact-candidate list 53 by the contact-candidate-list preparingunit 16 and the contact-force reference table 54 by thecontact-force-reference-table preparing unit 17 can be prepared throughthe following algorithm. “Contact-candidate-pair condition==true” in thefollowing algorithm means that Conditions a), b), and c) of thecontact-candidate list 53 are satisfied.

for(Search neighboring cell){   if(Contact-candidate-paircondition===true){     box_id = box_id+1     Lnum = s_jgi[i−1]+box_id    list_i[Lnum] = i     list_j[Lnum]= j     for(k = 0;k<n_jli[j];k++){      if(nei_pn[j][k]==i){         Ref[j][k] = Lnum       }     }   } }

In this embodiment, the contact-force reference table 54 includes onlyparticles that have a diameter greater than that of the target particle.Therefore, for the same reason as the neighboring-particle table 52described above, the amount of memory used for the sequences prepared inadvance for the contact-force reference table 54 is inhibited.

The contact-force-reference-table preparing unit 17 sends the derivedcontact-force reference table 54 Ref[i][box_id] to the contact-forcecalculating unit 19.

The contact determining unit 18 refers to the contact-candidate list 53to determine whether paired particles are in contact with each other forevery particle pair and calculates the interparticle contact force Fijrand Fijt in the case where the paired particles are in contact with eachother.

Specifically, upon receiving from the contact-candidate-list preparingunit 16 the contact-candidate list 53 list_i[Lnum] and list_j [Lnum],the integral contact-candidate number s_jgi[i], and the contact forceSpring[Lnum] of the contact-candidate pairs, which were also pairs inthe previous step, and the particle information or receiving from thecontact-candidate-list-update determining unit 13 the particleinformation, the contact determining unit 18 scans the contact-candidatelist 53 to determine the contact between the particles i and j of thecontact-candidate pairs. Particles i and j are determined to be incontact with each other, for example, if the calculated interparticledistance between particles i and j is smaller than or equal to the sumr_(i)+r_(j) of the diameters of the particles i and j.

If it is determined that the particles i and j are in contact with eachother, the contact force components Fn and Ft are calculated usingExpressions 1 and 2 of the above-described Voigt model. The tangentialcomponent of the current contact force of the spring is calculated usingExpression 2, where K_(t)x_(t) equals to spring[ ]+ΔFt, spring[ ]represents the tangential components of the contact force of the springin the previous step, as described above, and ΔFt represents an increase(or decrease) in the tangential component of the contact force of thespring from the previous step to the current step. ΔFt is calculated bymultiplying the relative speed between the particles i and j by the timeΔt of one step and the elastic modulus Kt of the tangential component ofthe spring. The elastic modulus Kn of the normal component of spring 72,the elastic modulus Kt of the tangential component of the spring 72 andthe coefficient of restitution e between the particles are values setprior to the simulation. The viscous damping coefficients ηn and ηt ofthe dashpot 73 are derived from Expression 3 using the given Kn, Kt, ande and masses mi and mj included in the particle information. Therelative interparticle displacement components xn and xt and therelative velocity components vn and vt are calculated on the basis ofthe particle coordinates (position information) and translationalvelocity and rotational velocity (velocity information) included in theparticle information. Fn and Ft are calculated from Expressions 1 and 2using the variables and parameters acquired as described above. Thecalculated Fn and Ft are used to calculate the translational componentFn+Ft and rotational component Ri×Fr of the contact force. Thesecalculated components are stored in sequences Fijt[Lnum] and Fijr[Lnum](sequences contained in the above-described contact-force boxes) havingthe contact-candidate-pair list numbers (contact-candidate numbers) Lnumas elements.

If the particles are not in contact, the contact force equals zero andthus zero is stored in the sequences Fijt[Lnum] and Fijr[Lnum](sequences contained in the above-described contact-force boxes) havingthe contact-candidate-pair list numbers Lnum as elements.

This contact force is calculated by the contact determining unit 18through the following algorithm.

for(Lnum){   i = list_i[Lnum]   j = list_j[Lnum]   if(Contactdetermination==true){     Calculation of contact force components (Fn,Ft)     Fijt[Lnum] = Fn + Ft     Fijr[Lnum] = Ri × Ft   } else {    Fijt[Lnum] = 0     Fijr[Lnum] = 0   } }

The contact determining unit 18 creates or updates the contact-forceboxes such that they are linked to the contact-candidate-pair-listnumbers Lnum every time the contact-candidate list 53 is created orupdated. Even if the contact-candidate list is not updated, thecontact-force boxes are created or updated using an existingcontact-candidate list upon the update of the particle information.

Contact determination and the contact force are calculated by threadingthe contact-candidate-pair-list numbers. The global memory of the GPUachieves optimal access efficiency when the access pattern is equivalentto the thread order. Since the particle pair numbers are random withrespect to the list numbers, the global memory in which the particleinformation is stored is accessed randomly, causing a significantdecrease in the memory-access speed. In contrast, a texture memoryhaving a cash feature can prevent a decrease in memory-access speed dueto random access. Therefore, the global memory region in which theparticle information is stored is set so that it can be referred to froma texture memory. Furthermore, since the particles are numbered inaccordance with the order of the cell number of the cells containing theparticles, the particles in a contact-candidate pair have close numbers.Moreover, the contact-candidate-pair list contains the particle numbersi in the ascending order. Therefore, by threading thecontact-candidate-pair-list numbers, the cash hit rate increases, andthe speed of processing routines for determining whether particlescontact each other and calculating the contact force is increased.

The contact determining unit 18 sends the information of the created orupdated contact-force boxes, together with the integralcontact-candidate numbers s_jgi[i] and the particle information, to thecontact-force calculating unit 19.

The contact-force calculating unit 19 calculates the translationalcomponent Ft_(all)[i] and rotational component Fr_(all)[i] of the totalcontact force acting on the particles i using the contact-forcereference table 54 and the integral contact-candidate numbers s_jgi.Upon receiving the contact-force reference table 54 from thecontact-force-reference-table preparing unit 17 and receiving theintegral contact-candidate numbers from the contact determining unit 18,the contact-force calculating unit 19 extracts the contact force betweeneach target particle i and a particle in contact with the particle ifrom the contact-force boxes using the contact-force reference table 54and the integral contact-candidate numbers s_jgi. Specifically, forparticles having a diameter greater than a target particle i, thecontact forces Fij[Lnum] corresponding to the box numbers are extractedfrom the contact-force boxes having box numbers (Lnum) (referenceinformation) stored in the contact-force reference table 54. For theparticles having a diameter smaller than that of the target particle i,the contact forces are extracted by specifying the storage positions ofthe contact-force tables using the integral contact-candidate numberss_jgi[ ]. More specifically, contact forces Fij[s_jgi[i−1]] toFij[s_jgi[i]−1] corresponding to the box numbers (contact-candidatenumbers) s_jgi[i−1] to s_jgi[i]−1 is extracted from the contact-forceboxes. The contact forces Fij[ ] stored in the contact-force boxes areeach stored as group of three one-dimensional sequences of thetranslational component Fijt[Lnum], the rotational component Fijr[Lnum],and the β[Lnum] used in Expression 13, which is described below.

The extraction of contact forces from contact force boxes will bedescribed in detail below with reference to FIG. 10. FIG. 10 is anexample in which the target particle i is numbered 4 (i=4), as in FIGS.6 and 9. The contact-force reference table 54 contains referenceinformation on the contact forces between the target particle (i=4) andparticles (i=7, 8) having a diameter greater than the target particle(i=4). The box number (contact-candidate-list number) Lnum of thecontact-force box containing the contact force between the targetparticle i=4 and the particle i=7 is 4, and the box number Lnum of thecontact-force box containing the contact force between the targetparticle i=4 and the particle i=8 is 7. The contact-force calculatingunit 19 uses these box numbers in the contact-force reference table 54to extract the contact force Fij[4] between the target particle i=4 andthe particle i=7 and the contact force Fij[7] between the targetparticle i 4 and the particle i=8 from the contact-force boxes.

As illustrated in FIG. 10, the integral contact-candidate numbers_jgi[i] of the target particle i=4, which is represented as s_jgi[4],is equal to 3. The integral contact-candidate number s_jgi[i−1] of theparticle having the particle number prior to the particle number of thetarget particle (which is i=3 in FIG. 10) is represented as s_jgi[3] andis equal to 1. Accordingly, for particles having a diameter smaller thanthe target particle i=4 (which are j=5, 6 in this example), the contactforces Fij[1] and Fij[2] are extracted from the contact-force boxesnumbered s_jgi[i−1] to s_jgi[i]−1, i.e., boxes numbered 1 and 2.

Then, the translational components FUN and the rotational componentsFr_(all)[i] of the all contact forces acting on every particle i arecalculated using the contact forces Fij[ ] (=Fijt[ ], Fijr[ ], β[ ])extracted from the contact-force boxes.

All contact forces (translational components) Ft_(all)[i] acting on aparticle i is obtained by summating the contact forces corresponding tothe contact-candidate-list numbers s_jgi[i−1] to s_jgi[i]−1 withoutchanging their signs. In contrast, the contact forces obtained byreferring to the contact-force reference tables 54 (Ref) numbered 0 ton_jli−1 are added together after reversing the signs. This isrepresented by the following expression.

$\begin{matrix}{{{Ft}_{all}\lbrack i\rbrack} = {{\sum\limits_{{box\_ id} = {{s\_ jgi}{\lbrack{i - 1}\rbrack}}}^{{{s\_ jgi}{\lbrack i\rbrack}} - 1}{{Fij}_{t}\lbrack{box\_ id}\rbrack}} - {\sum\limits_{{box\_ id} = 0}^{{{n\_ jli}{\lbrack i\rbrack}} - 1}{{Fij}_{t}\left\lbrack {{{Ref}\lbrack i\rbrack}\lbrack{box\_ id}\rbrack} \right\rbrack}}}} & (12)\end{matrix}$

All contact forces (rotary component) Fr_(all)[i] acting on the particlei are obtained through a similar addition process. To calculateFr_(all)[i], the contact forces corresponding to thecontact-candidate-list numbers s_jgi[i−1] to s_jgi[i]−1 are directlyreferred to, whereas the contact forces referred to in the contact-forcereference tables 54 (Ref) numbered 0 to n_jli−1 are summated after thesigns are reversed and β is multiplied. Fr_(all)[i] is determined by thefollowing expression, where β represents a value acquired by dividingthe distance from the center of a particle j to the contact point withthe distance from the center of a particle i to the contact point.

$\begin{matrix}{{{Fr}_{all}\lbrack i\rbrack} = {{\sum\limits_{{box\_ id} = {{s\_ jgi}{\lbrack{i - 1}\rbrack}}}^{{{s\_ jgi}{\lbrack i\rbrack}} - 1}{{Fij}_{r}\lbrack{box\_ id}\rbrack}} + {\sum\limits_{{box\_ id} = 0}^{{{n\_ jli}{\lbrack i\rbrack}} - 1}{{{Fij}_{r}\left\lbrack {{{Ref}\lbrack i\rbrack}\lbrack{box\_ id}\rbrack} \right\rbrack} \times {\beta \left\lbrack {{{Ref}\lbrack i\rbrack}\lbrack{box\_ id}\rbrack} \right\rbrack}}}}} & (13)\end{matrix}$

The contact forces corresponding to s_jgi[i−1] to s_jgi[i]−1 are contactforces from particles j to particles i, whereas the contact forcesreferred to in the contact-force reference tables 54 (Ref) numbered 0 ton_jli−1 are contact forces from the particles i to the particles j.Therefore, the translational components of the contact forces in thecontact-force reference tables 54 (Ref) numbered 0 to n_jli−1 aresummated after their signs are reversed, and the rotational componentsof the contact forces are summated after multiplying them by β[box id].By directly referring to the contact forces from the particles j to theparticles i in the sequences Fijt[ ] and Fijr[ ] containing the contactforces using s_jgi[i] but without using the contact-force referencetables 54 Ref[i][box id], the contact forces from numerous smallparticles to large particles, which can make contact-candidate pairswith the small particles, can be summated without an increase in memorysize used for the contact-force reference tables 54.

Specifically, through calculation by a computer, Ft_(all) and Fr_(all)are determined by threading particles i and running a loop of box_id.Since the memory randomly accesses Fij, a texture memory is used. Asdescribed above, since a texture memory has a high cash hit rate, theprocessing speed is greater than directly accessing the global memory.Continuous access to the contact-force reference tables 54 is possiblethrough array declaration to set the element size of B in the sequenceRef[A][B] is a factor of 16.

The contact-force calculating unit 19 sends, together with the particleinformation, the translational component Ft_(all)[i] and the rotationalcomponent Fr_(all)[i] of the total contact force acting upon eachparticle i, which are calculated as described above, to theparticle-information updating unit 20.

Upon receiving the translational components Ft_(all)[i] and therotational components Fr_(all)[i] of all contact forces acting on everyparticle i and the particle information from the contact-forcecalculating unit 19, the particle-information updating unit 20calculates the coordinates of the particle (position information) andthe translational and rotational velocities (velocity information) onthe basis of the positions and velocities of the particles and contactforces received from the contact-force calculating unit 19 and updatesthe particle information stored in the particle-information retainingunit 11 using the calculated coordinate, translational velocity, androtational velocity. Specifically, by solving a motion equation(Expression 4) of translational movement using the Leap-Flog method, thetranslational velocity and the particle coordinates at the currentmoment are determined. The motion equation (Expression 5) for rotationalmovement is also solved using the Leap-Flog method, and the acquiredrotational velocity is set as a predicted value. By calculating therolling friction corresponding to the predicted rotational velocity andsolving the motion equation of rotational movement in consideration ofthe rotating force and the rolling friction due to the contact force,the rotational velocity at the current moment is determined.

The particle-information updating unit 20 uses the calculated positioncoordinates, the translational velocity, and the rotational velocity ofeach particle to update the coordinates, the translational velocity, andthe rotational velocity in the particle information and sends theupdated particle information to the particle-information retaining unit11 to update the particle information in the particle-informationretaining unit 11.

With reference to FIG. 11, the processing by the particle simulator 10and the particle simulation method according to this embodiment will bedescribed below. FIG. 11 is a flow chart illustrating the processing bythe particle simulator 10 according to this embodiment. Theparticle-information acquiring unit 12 acquires particle information,which is input to and stored in the particle-information retaining unit11, from the particle-information retaining unit 11 (S101: Input step).The contact-candidate-list-update determining unit 13 determines whetherthe contact-candidate list is to be updated (S102). For example, thecontact-candidate list is updated if a particle of which the integraltraveling distance is greater than rα exists, where r represents theparticle radius and α represents a parameter representing the updatingfrequency. For example, α is set to 0.2. If updating is required, thenthe process proceeds to Step S103, else the process proceeds to StepS108.

Subsequently, the particle numbers are changed by the particle-numberchanging unit 14 (S103: Assigning step). The particle-number changingunit 14 renumbers the particles in accordance with the order of the cellnumbers of the cells containing the particles and renumbers theparticles in the same cell in accordance with the order of particlediameter. Then, the neighboring-particle-table preparing unit 15prepares the neighboring-particle table 52 including information ofneighboring particles of every particle in the work area 51 (S104).Particles that satisfy Conditions (a), (b), and (c) are stored in theneighboring-particle tables 52 as neighboring particles.

Subsequently, the contact-candidate-list preparing unit 16 calculatesthe integral contact-candidate numbers s_jgi[i], which are the prefixedsums of particle numbers n_jgi[i] of particles satisfying Condition (a)but not satisfying Conditions (b) and (c) of the neighboring-particletables 52 (S105: Integral contact-candidate-number calculating step).Then, the contact-candidate lists 53 including information of pairedparticles that may be in contact with each other among all particles inthe work area 51 are prepared (S106: Contact-candidate selecting step).Particles that satisfy Conditions a), b), and c) are stored in thecontact-candidate lists 53. The contact-force-reference-table preparingunit 17 prepares the contact-force reference tables 54 (S107:Contact-force reference table preparing step)

Subsequently, the contact determining unit 18 refers to thecontact-candidate lists 53 to determine whether two particles are incontact with each other (S108). If the particles are in contact witheach other, the interparticle contact force components Fijr and Fijt arecalculated using Expressions 1 and 2 and are stored in sequences ofcontact-force boxes (S109: Contact-force calculating step). If theparticles are not in contact with each other, the contact forcecomponents Fijr and Fijt, which equal to zero, are stored in thesequences of the contact-force boxes.

The contact-force calculating unit 19 refers to the contact-forcereference tables 54 and the integral contact-candidate numbers s_jgi[i]to select other particles in contact with target particles and obtainsthe contact forces from corresponding sequences Fijr and Fijt in thecontact-force boxes. The sum of the contact forces is calculated usingExpressions 12 and 13, and the contact force components Fr and Ft actingon every particle (S110: Contact-force-sum calculating step) arecalculated. The particle-information updating unit 20 calculates thecoordinates of every particle (position information) and thetranslational and rotational velocities (velocity information) on thebasis of the contact force components calculated in Step S109 and thecurrent particle information, and the particle information stored in theparticle-information retaining unit 11 is updated using the calculatedcoordinates, translational velocity, and rotational velocity (S111:Particle-information updating step).

Then, it is determined whether the particle simulation is completed(S112). For example, when a predetermined amount of time has elapsed, itis determined that the particle simulation is completed. When theparticle simulation is not yet completed, the process returns to StepS101 and proceeds to the next step. When the particle simulation iscompleted, the process ends.

The amount of memory used by the particle simulator 10 according to thisembodiment will be described below with reference to FIG. 12. FIG. 12illustrates the amounts of memory used in this embodiment (solid line)and comparative example (dotted line) at a variable ratio of the maximumparticle diameter to the minimum particle diameter of the particlemodel. The particle model includes 1,000,000 particles and is atwo-component system containing particles having a maximum particlediameter (maximum particles) and particles having a minimum particlediameter (minimum particles). The comparative example differs from theembodiment in that information on particles having a diameter smallerthan that of the target particle i are stored in theneighboring-particle table 52 and the contact-force reference table 54.In FIG. 12, the horizontal axis represents the ratio of the maximumparticle diameter to the minimum particle diameter, and the verticalaxis represents the memory size used (MB).

In the case of the comparative example, the memory size required for theneighboring-particle tables 52 nei_pn[i][box_id] and the contact-forcereference tables 54 Ref[i][box_id] is equal to the product of the numberof particles, the number of maximum contact-candidate particles, and 4bytes, if 4 bytes are needed for each data item. The number of maximumcontact-candidate particles is the maximum number of neighboringparticles stored in the neighboring-particle table 52 and thecontact-force reference table 54. In the case of the minimum particlesfilling the surface of a maximum particle, the number of maximumcontact-candidate particles is represented by the following expression:

4π(r _(max) +r _(min))² /πr _(min) ²  (14)

where r_(max) represents the diameter of a maximum particle, and r_(min)represents the diameter of a minimum particle. As illustrated in FIG.12, in the comparative example, as the maximum particle diameterincreases with respect to the minimum particle diameter, i.e., as theratio of the maximum particle diameter to the minimum particle diameterincreases, the number of maximum contact-candidate particles increases,causing an increase in the amount of memory used. For example, at aparticle diameter ratio of 1:10 (maximum particle diameter/minimumparticle diameter=10), the number of contact-candidate particles perparticle increases to 500 particles at maximum and requires 2 GB ofmemory.

For this embodiment, since particles having a diameter smaller than thatof the target particle are not selected as neighboring particles and arenot stored in the neighboring-particle tables 52 and the contact-forcereference tables 54, Expression 14 for calculating the number of maximumcontact-candidate particles does not hold. The number of maximumcontact-candidate particles in this embodiment does not depend on themaximum and minimum particle diameters. In this embodiment, the numberof maximum contact-candidate particles depends only on the parameter αincluded in dis_lim=(r_(i)+r_(j))(1.0+α), which is a threshold of theinterparticle distance in above-described Condition (a) of neighboringparticles but does not depend on the particle diameter distribution andis a constant. The amount of memory used is calculated for a case inwhich α is set to 0.2, and the number of maximum contact-candidateparticles is set to 16, which is a sufficient number satisfying anyparticle-diameter distribution condition. As a result, as illustrated inFIG. 12, the memory size used in this embodiment is a constant value andis equal to the product of the number of particles, the number ofmaximum contact-candidate particles, and 4 bytes, which equal 64 MB(1000000×16×4=64).

According to the particle simulator 10 and the particle simulationmethod according to this embodiment, for a particle neighboring aparticle i and having a diameter smaller than that of the particle i,among all particles in the work area 51, only the integralcontact-candidate numbers s_jgi[i], which represent the numbers ofparticles, is stored. The contact determining unit 18 uses the integralcontact-candidate number s_jgi[i] to specify a storage position in acontact-force table to extract the contact force between the particle iand another particle. In this way, the contact-force reference tables 54do not have to contain information associated with particles having adiameter smaller than that of particles i, among particles neighboringthe particles i, and merely information on particles having a diametergreater than that of the particles i needs to be stored. For simulationon a particle group having a particle diameter distribution, ifparticles having a diameter smaller than that of the particle areconsidered, the number of particles that may simultaneously contact theparticle will increases. Thus, a large memory size must be provided inadvance for the contact-force reference table 54. In contrast, if theabove-described configuration is used, merely the contact between theparticle and particles having a diameter greater than that of theparticle need to be considered. Thus, the number of particles for whichcontact determination need to be performed is decreased, and the amountof memory required for preparing the contact-force reference tables 54can be decreased. As a result, the memory size used can be inhibited. Byinhibiting the memory size used, particle simulation can be performedusing a shared-memory parallel computer, such as a GPU.

The particle-number changing unit 14 specifies all particles in the workarea 51 in the order corresponding to the cell numbers, assigns particlenumbers to the specified particles contained in the same cell in theorder corresponding to the particle diameter, and sorts the particles inthe order of the cell numbers and, for particles in the same cell, inthe order of the particle diameter. Accordingly, during preparation ofthe contact-force reference tables, particles having a small diametercan be easily recognized on the basis of the particle numbers and do nothave to be selected for entry to the contact-force reference tables. Asa result, the calculation efficiency can be improved for preparation ofthe contact-force reference table.

An embodiment of the present invention has been described above. Thepresent invention, however, is not limited to this embodiment. In theabove-described embodiment, the particle-number changing unit 14 sortsthe particles in two steps: 1) the particles are arranged in accordancewith the order of the cell numbers containing the particles; and 2) theparticles in the same cell are rearranged in the order of particlediameter. Instead of performing the two steps, only the first step maybe performed.

In the above-described embodiment, a GPU is used for the presentinvention. Instead, shared-memory hardware (for example, a vectorprocessor) that performs single instruction multiple data (SIMD)parallel processing in the same manner as the GPU, i.e., a shared-memoryparallel computer, may be used.

1. A particle simulator simulating the behavior of particles in aparticle group having a diameter distribution in a work area throughcalculation of positions and velocities of the particles based on valuesof contact forces generated between the particles and other particles inthe particle group, the particle simulator comprising: input means forinputting particle information including position information andvelocity information of the particles; particle-information storingmeans for storing the particle information; assigning means forassigning specific information items to the particles to identify theparticles in an order in accordance with the position information;contact-candidate selecting means for selecting target particles fromparticle group in the ascending order of the specific information itemsand selecting non-target particles paired with the target particles ascontact-candidate pairs, the non-target particles presumably being incontact with the target particles and having a diameter smaller than thediameter of the target particles; contact-force calculating means forcalculating the values of contact forces generated between the twoparticles in the contact-candidate pairs selected by thecontact-candidate selecting means based on the position information andthe velocity information of the paired particles and storing informationitems of the calculated values of contact forces in an ordercorresponding to the order of the contact-candidate pairs incontact-force tables; contact-force-reference-table preparing means forpreparing contact-force-reference tables including reference informationused to refer to the contact-force tables to retrieve the values ofcontact forces generated between the particles and neighboring particlesneighboring the particles and having a diameter greater than thediameter of the particles; integral-contact-candidate-number calculatingmeans for calculating integral contact-candidate numbers being integralvalues of contact-candidate numbers indicating the numbers ofneighboring particles neighboring the particles having a diametersmaller than the diameter of the particles; contact-force-sum computingmeans for extracting the values of contact forces generated between theparticles and other particles having a diameter greater than thediameter of the particles by using the reference information in thecontact-force-reference tables, extracting the values of contact forcesgenerated between the particles and other particles having a diametersmaller than the diameter of the particles by specifying the storagepositions of the values of contact forces in the contact-force tablesusing the integral contact-candidate numbers, and computing the sum ofcontact forces for each particles using the extracted values of contactforce; calculating means for calculating new position information andvelocity information of the particles based on the values of contactforces of the particles calculated by the contact-force-sum computingmeans; and particle-information updating means for updating the particleinformation stored in the particle-information storing means using thenew position information and velocity information calculated by thecalculating means.
 2. The particle simulator according to claim 1,wherein the work area is segmented into a plurality of cells, each cellhaving a cell number, the particle information includes the cell numbersof the cells specified by the position information and containing theparticles, and the assigning means identifies the particles in an ordercorresponding to the cell numbers and assigns specific information itemsspecified in accordance with the order corresponding to the particlediameter to particles in the same cell.
 3. The particle simulatoraccording claim 2, wherein the contact-candidate selecting means selectsthe target particles from the particle group in the ascending order ofthe specific information items and prepares contact-candidate lists ofthe target particles including, for all of the particles in the workarea, contact-candidate-pair information including pairs of the targetparticles and the non-target particles presumably being in contact withthe target particles and having a diameter smaller than the diameter ofthe target particles and contact-candidate numbers assigned to the pairsin the ascending order; the contact-force calculating means calculatesthe values of contact forces generated between the paired particlescorresponding to the contact-candidate-pair information included in thecontact-candidate lists are calculated based on the position informationand the velocity information of the paired particles and storescontact-force information linking the calculated values of contactforces and the contact-candidate number in contact-force tables; thecontact-force-reference-table preparing means prepares thecontact-force-reference tables linking the contact-candidate numbersassociated with the particles and other particles with the particles andthe other particles as reference information used to refer to thecontact-force tables to retrieve the values of contact forces generatedbetween the particles and the other particles presumably being incontact with the particles and having a diameter greater than thediameter of the particles; the integral-contact-candidate-numbercalculating means calculates integral contact-candidate numbers s_jgi[i](where i represents specific information items of the particles) beingintegral values integrating, in the order of the specific informationitems, contact-candidate numbers indicating the number of otherparticles presumably being in contact with the particles and having adiameter smaller than the diameter of the particles; and thecontact-force-sum computing means extracts the values of contact forcesof the particles from the contact-force tables using the referenceinformation of the contact-force-reference tables, extracts the valuesof contact forces corresponding to the contact-candidate numberss_jgi[i−1] to s_jgi[i]−1 from the contact-force tables using theintegral contact-candidate numbers s_jgi[i], and calculates the sum ofthe contact forces of the particles using the extracted information ofthe contact forces.
 4. A method of simulating particles in a particlesimulator, the method simulating the behavior of particles in a particlegroup having a diameter distribution in a work area by calculating thepositions and velocities of the particles based on values of contactforces generated between the particles and other particles in theparticle group, the method comprising, the steps of: inputting particleinformation including position information and velocity information ofthe particles and storing the particle information inparticle-information storing means; assigning specific information itemsto the particles to identify the particles in an order in accordancewith the position information; selecting target particles from particlegroup in the ascending order of the specific information items andselecting non-target particles paired with the target particles ascontact-candidate pairs, the non-target particles presumably being incontact with the target particles and having a diameter smaller than thediameter of the target particles; calculating the values of contactforces generated between the two particles in the contact-candidatepairs selected in the selecting step based on the position informationand the velocity information of the paired particles and storinginformation items of the calculated values of contact forces in an ordercorresponding to the order of the contact-candidate pairs incontact-force tables; preparing contact-force-reference tables includingreference information used to refer to the contact-force tables toretrieve the values of contact forces generated between the particlesand neighboring particles neighboring the particles and having adiameter greater than the diameter of the particles; calculatingintegral contact-candidate numbers being integral values ofcontact-candidate numbers indicating the numbers of neighboringparticles neighboring the particles having a diameter smaller than thediameter of the particles; extracting the values of contact forcesgenerated between the particles and other particles having a diametergreater than the diameter of the particles by using the referenceinformation in the contact-force-reference tables, extracting the valuesof contact forces generated between the particles and other particleshaving a diameter smaller than the diameter of the particles byspecifying the storage positions of the values of contact forces in thecontact-force tables using the integral contact-candidate numbers, andcomputing the sum of contact forces for each particles using theextracted values of contact force; calculating new position informationand velocity information of the particles based on the values of contactforces of the particles calculated in the extracting step; and updatingthe particle information stored in the particle-information storingmeans using the new position information and velocity informationcalculated by the calculating step.
 5. The method of simulatingparticles according to claim 4, wherein the work area is segmented intoa plurality of cells, each cell having a cell number; the particleinformation includes the cell numbers of the cells specified by theposition information and containing the particles; and in the assigningstep, the particles in an order corresponding to the cell numbers areidentified and assign specific information items in accordance with theorder corresponding to the particle diameter to particles in the samecell.